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Abstract. A method for the calculation of the temperature dependence of dielectric 
and piezoelectric responses, based on the use of a first-principles effective Hamiltonian, 
is described. Results are presented for the ferroelectric perovskite PbTi0 3 . While the 
method includes only the soft-mode contributions to the responses, it is argued to give 
a good description of the divergences or near-divergences of the response functions 
near the cubic-tetragonal transition. The expression of the response functions in terms 
of correlation functions is used to provide a real-space interpretation of the responses 
which clearly distinguishes between PbTi03 and the related materials BaTiC>3 and 
KNb0 3 . 



INTRODUCTION 

Systems with large dielectric and piezoelectric responses are of interest both from 
a fundamental and a technological point of view [1] . Understanding the microsopic 
origin of a high sensitivity of the polarization of a ferroelectric oxide to applied 
field, in the case of the dielectric constant, or stress- induced changes in strain, in 
the case of piezoelectricity, can lead to the optimization of these properties through 
the appropriate choice of materials. Empirically, large responses are observed in 
stoichiometric compounds such as PbTi03 and BaTi03 for temperatures near the 
ferroelectric phase transition [2-5]. In solid solutions, such as PZT, enhanced re- 
sponses are observed over a wide temperature range [6]. 

In this paper, we discuss the modeling of the piezoelectric and dielectric re- 
sponse of ferroelectric perovskites, specifically PbTi0 3 and BaTi0 3 , though the 
use of first-principles effective Hamiltonians previously constructed in Refs. [7-10]. 
The calculations give reasonable agreement with experiment, considering that this 
modeling gives only the contribution of the soft modes to the responses, and not 
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other contributions such as thermal expansion and the response of other polar 
modes. However, for temperatures near the ferroelectric phase transition, the soft 
modes are expected to dominate and thus this approach should capture the essen- 
tial physics of the large observed responses. In addition, it should correctly describe 
the trends between different materials. Finally, with this microscopic approach, we 
are able directly to relate the calculated response functions to the characteristic 
correlations of local unit cell polarizations, which should yield further insight into 
the nature of the large responses of interest. 



FIRST-PRINCIPLES EFFECTIVE HAMILTONIANS 

The construction of first-principles effective Hamiltonians for PbTi0 3 and 
BaTiOs has previously been described in detail [8-10]. Briefly, the model con- 
sists of one vector per five-atom unit cell which represents the local polarization 
associated with the ferroelectric distortion. The potential 7~t e ff({£,i}, ^ap) for these 
vector degrees of freedom is expanded for the high-symmetry cubic perovskite ref- 
erence structure, including local anharmonic terms, quadratic intersite interactions 
which are assumed to be dipolar beyond third neighbors, and lowest order coupling 
to homogeneous strain. 

The extension of this model to include the effects of electric field is accomplished 
by writing the dependence of the polarization on the model degrees of freedom to 
lowest order: 

P({£i},e afj ) = '£Tea< ) £ i (1) 

i 

and adding to the potential the coupling term — P ■ E where E is the macroscopic 
electric field [11]. The dielectric response is then simply obtained from d ^° > , 

while the piezoelectric response can be expressed as d< ^ p> , where the brackets are 
used to denote the thermal expectation value. The effects of macroscopic stress 
oij can be included by adding to the potential the coupling term — cr^- • e^. With 
these couplings, the thermodynamic identity = jj^- can be obtained (the well- 
known equality of the direct and converse piezoelectric effects [11]). In the rest of 
this paper, we work at zero external stress and evaluate the piezoelectric response 

• d<e a g> 
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calculations of the temperature dependent properties of the system are carried 
out using classical Monte Carlo calculations, as previously described [8,9]. The cal- 
culations of the response functions presented here involved 7x7x7 simulation cells, 
10,000 Monte Carlo sweeps (MCS) for thermalization, 200,000 MCS for the com- 
putation of thermal expectation values, and several runs at each temperature with 
different random number seeds to estimate statistical error. The calculations of the 
real-space correlation functions presented here involved 10x10x10 simulation cells, 
10,000 Monte Carlo sweeps (MCS) for thermalization, and 100,000 MCS for the 



computation of thermal expectation values. In order to compute the temperature- 
dependence of the order parameter and thus identify the cubic-tetragonal transi- 
tion, a small symmetry-breaking field E z of magnitude ^- kV/cm was applied. The 
response function results include this nonzero field. 



CORRELATION FUNCTION EXPRESSIONS FOR 
RESPONSE FUNCTIONS 

In the Monte Carlo simulations described above, the response functions are com- 
puted using the fact that they can be expressed as correlation functions. This 
allows the appropriate derivatives to be calculated in one Monte Carlo run. 

Specifically, we consider the dielectric tensor: 

e«/3 = Coo + 47TXa/3 

where 

_ d < fr^Pg > 

Here, fl is the volume of the unit cell and P is the polarization per unit cell, 
specifically 

, Rta y - Si^MUj 4)( ^ St Uexp(-(3(H eff - Z*ea &, ■ E)) 
f{de a p}{Uj d£j}exp(-(3(H eff - Z*ea J2i 6a • E)) 



or 
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'N(3' dE a 
where we have defined 

Z((3,E) = f{de aP }{l[di j }ex P (-p(H eff -Z*eaoJ2^a-E)). 

J 3 i 

Using the approximation d<n 9 ^ a> = 9< d ^ (see footnote 7 in Ref. [12]) and 
differentiating Equation 2 with respect to E a , we readily find the correlation func- 
tion expression: 

x a , = ^pl(< £^£^ >-n<u >< & >) (3) 

» 3 

where 



The piezoelectric tensor can be expressed in terms of correlation functions in an 
completely analogous way: 

_ d < e a/3 > 
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where 



S{de a 't3>}{Ylj d£j}e af3 exp(-p(H eff - Z*ea £j & • E)) 

< e ad >= ~ =: ~ • (4) 

J{de a/ ^}{Ujd^}exp(-P(H eff - Z*ea J2 l & • E)) 

Differentiating Equation 4 with respect to E 7 , we find the correlation function 
expression: 

djap = (3(Z*ea )(< e a/3 J2^n > -N < e a/3 >< £ 7 >). (5) 



When the homogeneous strain appears only up to quadratic order in the effective 
Hamiltonian, as is the case for the effective Hamiltonians available for PbTi0 3 [8,9], 
BaTiOa [10] and KNbC>3 [13], the piezoelectric tensor can be reexpressed in terms of 
correlations of local polar distortions. Because the integral over strain is Gaussian, 
completely equivalent expressions for the thermal expectation values involving e a p 
can be obtained by making the following substitutions (and cyclic permutations 
thereof): 

- • lf ( " + f) El£i 2 + ^ 91 n , E(2& -&-&)) (6) 



N\C ll + 2C l2 )^^ 3(C n -C 12 ) 



e 



1 92 

xy ^ N C 44 



For the example of (i 33 , this yields 



+ ~,^ - ^ x £(2(< iUjz > -<zl >< ijz >) 

"(< fj^* > - < £y >< tjz >) 
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While Equation 5 proves more convenient for the actual computation of rf 7Q/ 3, it 
is Equation 8 that naturally leads to a microscopic interpretation of the soft-mode 
response, as we will discuss further below. 



e zz AND d 33 FOR PbTiQ 3 



The results of the Monte Carlo calculations for e zz are shown in Figure 1. At 
room temperature, the calculated value is 66, which should be compared to the 
single-crystal measurement e zz = 80 [14]. The agreement is surprisingly good, 
considering that the calculation includes only the soft-mode contribution. The 
room temperature value is already considerably larger than the zero-temperature 
calculated value of 42. The overall temperature dependence is very similar to that 
observed experimentally [2,14], with a slow increase with temperature below the 
transition, a near-divergence both above and below T c which is cut off by the first- 
order character of the cubic-tetragonal transition, and an enhanced value of the 
dielectric response in the high-temperature cubic phase. 

The results of the Monte Carlo calculations for ^33 are shown in Figure 2. At room 
temperature, the calculated value is 76, which should be compared to the single- 
crystal measurement rf 33 = 83.7 [14]. As for the dielectric response, the agreement 
is surprisingly good, considering that the calculation includes only the soft-mode 
contribution. The room temperature value is already considerably larger than 
the zero-temperature calculated value of 48. The overall temperature dependence 
is very similar to that observed experimentally [2,4], with a slow increase with 
temperature below the transition, a near-divergence below T c which is cut off by 
the first-order character of the cubic-tetragonal transition, and a drop to zero in 
the high-temperature cubic phase, as required by symmetry. 

These results can be compared with the closely analogous calculation of the 
piezoelectric coefficients of BaTiOs (Figure 1 in [15]), noting that for this system 
as well, remarkably good quantitative agreement with the experimentally measured 
values is obtained. 
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FIGURE 1. Dielectric response e zz as a function of temperature, in K. 



LANDAU THEORY FOR THE DIVERGENT 

RESPONSES 

The near-divergences observed in the results described above can be already 
understood within a simple Landau theory which describes coupling of a scalar 
polar mode £ with a one-component strain e and electric field E: 



F{£, e; T, E) = A£ 2 + 5£ 4 - E£ + ce 2 + ge? 



(9) 



where the temperature dependence enters through A = A (T — T c ). If we first 
eliminate e, which is accomplished by the substitution e — * ^£ 2 , the following 
effective fourth-order free energy results: 



F(£,e;T,E) = At? + (B - ^)£ 4 - E£ 



(10) 



Minimizing with respect to £ yields the following expressions for the dielectric 
susceptibility 
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and piezoelectric coefficient 
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FIGURE 2. Piezoelectric response ^33 as a function of temperature, in K. 



As observed in the full simulations, d diverges as T c is approached from below, 
while x diverges on both the low and high temperature sides, and is larger above 
T c . In PbTiOs, these divergences are cut off due to the first-order character of the 
transition, so that, for example, a Curie- Weiss fit to x above T c yields a value of 
T about 80 K below T c . 



REAL-SPACE DECOMPOSITION OF THE 
CORRELATION FUNCTIONS 

The correlation function expressions for the response functions allow the interpre- 
tation of the responses in terms of the characteristic intersite correlations between 
the local polar distortions. Some of the correlation functions of interest have been 
previously calculated in effective Hamiltonian studies of KNb0 3 [13] and BaTi0 3 
[16], allowing us to connect to and build upon the insights derived in those inves- 
tigations to develop a microscopic understanding of the dielectric and piezoelectric 
responses. 

For the dielectric response, the first term inside the parentheses of Equation 3 
can be reorganized: 

jr £(< tiatif) > + E < > + E < &a&8 > +"0 

i j,nn j,nnn 

If correlations become infinite ranged, as at a second order transition, this will scale 
like N, leading to a divergence in the thermodynamic limit. 

For the piezoelectric tensor, we can similarly reorganize the correlation functions 
which appear in Equation 8, for example: 

i j,nn j,nnn 

The real-space correlation functions appearing in the expressions above are com- 
puted using fast fourier transforms of & and Ciadp- For the dielectric susceptibility, 
for example, the real space expressions can be written in terms of Xap(o) as follows 
(taking Rj = Ri + a): 

E < &a&/3 >= ^ E E < > exp(-iq- Ri-iq'- (Ri + d)) 

i i q,q ' 

= E < Ca(q)^(q)* > exp(iq- d) 
It can be checked that this gives the correct relation to Xa/3 when < £ >= 0: 
Xap = ^ Z ^ E E < tictiP >= E E Xaf3(q)exp(iq- d) = Xaf3(q = 0) 

d i d Q 



In Figure 3, we show the real-space decomposition in the xz plane of the high- 
T dielectric susceptibility correlation function < £,i Z £ i+ j z > f° r PbTi0 3 at two 
temperatures above T c . The correlation is rather isotropic, and the range of the 
correlations increases as expected as the temperature decreases towards T c . For 
comparison, the same correlation function for BaTiC>3 at two temperatures above 
the cubic-tetragonal transition temperature is also shown in Figure 3. The highly 
anisotropic "chain-like" correlations discussed in previous work [13,16] are clearly 
evident, and the correlation length transverse to the chains increases as the tem- 
perature decreases towards T c . Similar behavior is expected for KNb0 3 . The 
qualitative difference in the nature of the correlations can be directly attributed to 
the difference in the dispersion relation of the soft mode along R. The implications 
of the different character of these correlations for the dielectric behavior of PbTiC>3 
and BaTiC>3 are currently under investigation. 

As can be seen from Equation 8, the piezoelectric response depends on third- 
order two-site correlations of the form < > ~~ < Cf a >< CjP >■ These have 
not to our knowledge been considered in previous studies. Calculations of these 
real-space correlations for PbTi0 3 and BaTi0 3 above and below the transition 
temperature are currently in progress. 
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FIGURE 3. Real-space decomposition in the xz plane of the high-T dielectric susceptibility 
correlation function < £iz€ i+ j tZ > for (a) PbTi0 3 at T = 706 K, (b) PbTi0 3 at T = 690 K, (c) 
BaTi0 3 at T = 343 K, and (d) BaTi0 3 at T = 308 K. The numerical value in the box indicates 
the on-site correlation< £i Z £i Z >, while the numbers above and to the right of the box give the 
correlations for the corresponding values of d. All correlations have been normalized by dividing 
by the square of the ground state value of £ 2 (0.08 2 for PbTi0 3 and 0.03 2 for BaTi0 3 ). Only the 
largest correlations are shown (above 0.0078 for PbTi0 3 and above 0.056 for BaTi0 3 ). 



DISCUSSION AND CONCLUSIONS 



As already discussed, the use of an effective Hamiltonian for calculating the 
dielectric and piezoelectric response involves a number of approximations. The 
effects of thermal expansion and the neglect of other polar modes cannot be incor- 
porated without adding additional degrees of freedom to the model (unless one is 
willing to accept semi-empirical input). However, the low-order expansion around 
the cubic reference structure can be improved, which in particular will yield bet- 
ter model results for the properties of the PbTiOs tetragonal ground state [19]. 
This could involve both H e ff and P. In the case of the polarization P, the next 
order terms involve including strain dependence and £ dependence of the effective 
charges. First-principles calculations [17,18] suggest that the former corrections are 
rather small compared to the latter. Calculations for these refinements for PbTiOs 
are currently in progress. 

Since the effective Hamiltonian approach is limited to the description of soft-mode 
contributions to the responses, we focus on the near-divergent responses where these 
contributions strongly dominate. Fortunately, this regime is, in any case, the one 
of greatest interest for understanding the origin of and engineering large responses 
in a variety of systems characterized by proximity to a lattice instability. Further 
insight can be gained from the real-space decomposition of the divergent responses 
based on their expressions in terms of correlation functions. 

Finally, these results for pure systems suggest an interpretation for the origin 
of large responses in solid solutions. For such systems, there is a distribution of 
local environments, and for a given temperature some subset of lattice degrees 
of freedom will be marginally stable. This idea has been investigated using first 
principles results in Pbi^Ge^Te, and more details are presented in Refs. [20,21]. 
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